library("terra")
library("rstac")
options(timeout = 600)Skrypty geoprzetwarzania
Klasteryzacja danych
Pozyskanie danych
Niniejszą analizę wykonamy na podstawie danych z Sentinela-2. Jako źródło danych wykorzystamy katalog STAC i usługę Earth Search. Przeprowadzimy również proste filtrowanie dostępnych produktów definiując w zapytaniu: kolekcję (collections), zakres przestrzenny w układzie WGS 84 (bbox), interwał czasowy w standardzie RFC 3339 (datetime) oraz atrybut zachmurzenia sceny (eo:cloud_cover).
stac_source = stac("https://earth-search.aws.element84.com/v1")
stac_source |>
stac_search(
collections = "sentinel-2-c1-l2a",
bbox = c(22.5, 51.1, 22.6, 51.2),
datetime = "2023-03-01T00:00:00Z/2023-10-31T00:00:00Z") |>
ext_query(`eo:cloud_cover` < 10) |>
post_request() -> obrazyJako wynik zapytania zostały zwrócone metadane scen spełniające zadane warunki. Następnie możemy wybrać przykładową scenę z 21 września 2023 r. oraz ograniczyć liczbę kanałów do czterech podstawowych w rozdzielczości 10 m, tj. niebieskiego, zielonego, czerwonego i bliskiej podczerwieni w celu uproszczenia analizy. Finalnie, używając funkcji assets_url() zostaną zwrócone odpowiednie odnośniki do pobrania rastrów.
kanaly = c("blue", "green", "red", "nir")obrazy |>
items_filter(properties$`s2:tile_id` == "S2A_OPER_MSI_L2A_TL_2APS_20230921T151458_A043076_T34UEB_N05.09") |>
assets_select(asset_names = kanaly) |>
assets_url() -> sentinel
sentinel[1] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B02.tif"
[2] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B03.tif"
[3] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B08.tif"
[4] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B04.tif"Zauważ, że kanał 8 (bliska podczerwień) jest przed kanałem 4 (czerwony). Należy mieć na uwadze, że nieodpowiednia kolejność może spowodować problemy w dalszej części projektu.
W kolejnym kroku, używając funkcji dir.create() stworzymy nowy katalog na dysku, do którego zostaną pobrane zobrazowania. Oprócz tego, należy zdefiniować ścieżki i nazwy plików. W tym celu będą przydatne dwie funkcje – file.path() do stworzenia ścieżek do plików oraz basename() do wyodrębnienia nazwy pliku wraz z rozszerzeniem z URL.
dir.create("sentinel")
rastry = file.path("sentinel", basename(sentinel))Teraz możemy pobrać nasze dane używając funkcji download.file() w pętli.
for (i in seq_along(sentinel)) {
download.file(sentinel[i], rastry[i], mode = "wb")
}Przygotowanie danych
Po pobraniu danych musimy stworzyć listę plików (rastrów), które zamierzamy wczytać. W tym celu możemy wykorzystać funkcję list.files(), która jako argument przyjmuje ścieżkę do folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy wczytać (pattern = "\\.tif$") oraz zwrócić pełne ścieżki do plików (full.names = TRUE).
rastry = list.files("sentinel", pattern = "\\.tif$", full.names = TRUE)
rastry[1] "sentinel/B02.tif" "sentinel/B03.tif" "sentinel/B04.tif" "sentinel/B08.tif"Kiedy utworzyliśmy już listę plików, możemy je wczytać za pomocą funkcji rast(). Pobrane zobrazowania pokrywają duży obszar (około 12 000 km\(^2\)). Dla uproszczenia analizy zdefiniujmy mniejszy zakres przestrzenny do wczytania używając funkcji ext() oraz przekazując obiekt SpatExtent jako argument win w funkcji rast().
Zwróć uwagę, że do zdefiniowania zakresu przestrzennego podczas wyszukiwania danych użyliśmy układu WGS 84, natomiast teraz wymagany jest rzeczywisty układ rastra, tj. EPSG:32634. Można to sprawdzić w metadanych katalogu STAC (obiekt obrazy) lub używając funkcji describe() (wymaga ścieżki do pliku).
bbox = ext(510000, 540000, 5630000, 5650000)
r = rast(rastry, win = bbox)Możemy również zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.
names(r) = kanaly
rclass : SpatRaster
dimensions : 2000, 3000, 4 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
window : 510000, 540000, 5630000, 5650000 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 34N (EPSG:32634)
sources : B02.tif
B03.tif
B04.tif
B08.tif
names : blue, green, red, nir
min values : > 1, > 352, > 677, > 910
max values : 18560<, 17648<, 17056<, 16514< Następnie możemy sprawdzić podstawowe statystyki opisowe używając funkcji summary(). W przypadku dużych rastrów, zostanie automatycznie wykorzystana próba 100 tys. komórek.
summary(r)Warning: [summary] used a sample blue green red nir
Min. :-0.00080 Min. :0.00440 Min. :0.0027 Min. :0.0008
1st Qu.: 0.01510 1st Qu.:0.03260 1st Qu.:0.0211 1st Qu.:0.1951
Median : 0.03130 Median :0.06070 Median :0.0497 Median :0.2371
Mean : 0.03579 Mean :0.06165 Mean :0.0622 Mean :0.2491
3rd Qu.: 0.05170 3rd Qu.:0.08420 3rd Qu.:0.0958 3rd Qu.:0.2928
Max. : 0.49840 Max. :0.49280 Max. :0.5080 Max. :0.6909 Zasadniczo, wartości odbicia spektralnego mieszczą się w przedziale od 0 do 1, gdzie 0 oznacza brak odbicia (całe światło zostało pochłonięte), a 1 oznacza całkowite odbicie od powierzchni. W praktyce obiekty nie odbijają bądź pochłaniają stuprocentowo światła, niemniej sensory oraz procesy kalibracyjne nie są idealne, więc mogą pojawić się wartości odstające od tego przedziału, tak jak w tej sytuacji.
Można ten problem rozwiązać na dwa sposoby:
- Zastąpić te wartości brakiem danych (
NA). - Dociąć do minimalnej i maksymalnej wartości.
Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.
# sposób 1
r = clamp(r, lower = 0, upper = 1, values = FALSE)# sposób 2
r = clamp(r, lower = 0, upper = 1, values = TRUE)Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym przypadku zamiast funkcji plot() należy użyć funkcji plotRGB() oraz zdefiniować kolejność kanałów czerwonego, zielonego oraz niebieskiego. Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować rozciągnięcie kolorów używając argumentu stretch = "lin" lub stretch = "hist".
# plotRGB(r, r = 3, g = 2, b = 1)
plotRGB(r, r = 3, g = 2, b = 1, stretch = "lin")
Klasteryzacja
library("cluster") # klasteryzacja danychDane do modelowania muszą zostać przygotowane w odpowiedni sposób. Modele klasyfikacyjne najczęściej na etapie trenowania wymagają macierzy lub ramki danych (data frame). Dane rastrowe można przetworzyć do macierzy przy użyciu funkcji values().
mat = values(r)
nrow(mat) # wyświetla liczbę wierszy[1] 6000000Za pomocą interaktywnej funkcji View() możemy sprawdzić jak wygląda nasza macierz.
View(mat)W macierzy występują brakujące wartości. Zazwyczaj modele nie obsługują NA, więc musimy je usunąć. Służy do tego dedykowana funkcja na.omit().
mat_omit = na.omit(mat)
nrow(mat_omit)[1] 5999675Teraz przejdziemy do kolejnego etapu analizy, jakim jest wytrenowanie modelu. Istnieje wiele metod i modeli grupowania (patrz CRAN Task View), ale w tym przykładzie użyjemy prostego modelu grupowania metodą k-średnich. Ten model wymaga jedynie, aby podać z góry liczbę grup/klastrów (argument centers). Jest to algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby analiza była powtarzalna musimy ustawić ziarno losowości – set.seed().
set.seed(123) # ziarno losowości
mdl = kmeans(mat_omit, centers = 3)W wyniku powyższej operacji otrzymaliśmy m.in.:
- Obliczone średnie wartości grup dla poszczególnych kanałów (
mdl$centers). - Wektor ze sklasyfikowanymi wartościami macierzy (
mdl$cluster).
Wyświetlmy te obiekty:
mdl$centers blue green red nir
1 0.06072469 0.09172913 0.11470107 0.2175304
2 0.01790128 0.03401041 0.02827088 0.1998795
3 0.02791346 0.05896356 0.04148333 0.3397331head(mdl$cluster)[1] 2 2 2 2 2 2Oznacza to, że pierwszy wiersz (reprezentujący pojedyncze oczko siatki) należy do grupy 2, drugi do grupy 2, trzeci do grupy 2, itd.
Walidacja
Nieodłącznym elementem modelowania jest walidacja opracowanych modeli. Wyzwaniem jest wybór właściwej metody grupowania dla konkretnego zbioru danych i określenie odpowiedniej liczby grup. Należy pamiętać, że zwiększenie liczby klastrów zwiększa podobieństwo między obiektami w klastrze, ale przy ich większej liczbie interpretacja staje się trudniejsza.
Najczęstszym sposobem walidacji wyników grupowania jest użycie wewnętrznych metryk, takich jak wskaźnik Dunna, wskaźnik Daviesa-Bouldina lub wskaźnik sylwetki (silhouette index). Na potrzebny niniejszej analizy użyjemy tego ostatniego.
Indeks sylwetki ocenia zbieżność i separację klastrów na podstawie odległości między obiektami w tym samym klastrze i między obiektami w różnych klastrach. Wartości tego wskaźnika mieszczą się w zakresie od -1 do 1. Wartość bliska 1 wskazuje, że obiekt jest dobrze zgrupowany i znajduje się daleko od sąsiednich klastrów. Wartość bliska -1 sugeruje, że obiekt mógł zostać przypisany do niewłaściwego klastra. Wartość bliska 0 wskazuje, że obiekt znajduje się bardzo blisko granicy pomiędzy różnymi klastrami. Ogólnie rzecz biorąc, wyższa wartość tego wskaźnika wskazuje na lepsze wyniki grupowania. Więcej szczegółów można znaleźć w dokumentacji cluster::silhouette().
Spróbujmy teraz obliczyć wartości tego wskaźnika. Zasadniczo wymaga to obliczenie podobieństwa każdego obiektu do każdego obiektu, co w naszym przypadku jest zadaniem niemożliwym (nasz zbiór danych składa się z ponad 6 mln obiektów). Aby to wykonać, musimy wykorzystać mniejszą próbkę (załóżmy \(n=10000\)). W funkcji musimy określić dwa obiekty –- wektor z klastrami oraz macierz niepodobieństwa, którą można wcześniej obliczyć za pomocą funkcji dist().
set.seed(123)
# losowanie indeksów
idx = sample(1:nrow(mat_omit), size = 10000)
head(idx)[1] 2017743 3334670 4328362 3269750 2782437 5891675# obliczenie wskaźnika sylwetki
sil = silhouette(mdl$cluster[idx], dist(mat_omit[idx, ]))
summary(sil)Silhouette of 10000 units in 3 clusters from silhouette.default(x = mdl$cluster[idx], dist = dist(mat_omit[idx, ])) :
Cluster sizes and average silhouette widths:
3450 3423 3127
0.4483002 0.4448400 0.4106674
Individual silhouette widths:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.0601 0.3039 0.4940 0.4353 0.5926 0.6567 Średnia zbieżność klastrów wynosi 0,44. Nie jest to najlepszy wynik (powinniśmy spróbować zwiększyć liczbę klastrów lub użyć innej metody grupowania). Możemy również zaprezentować wyniki na wykresie.
kolory = rainbow(3) # wybierz 3 kolory z wbudowanej palety `rainbow`
plot(sil, border = NA, col = kolory, main = "Silhouette Index")
Interpretacja
Istotą grupowania jest utworzenie grupy podobnych obiektów, natomiast naszym zadaniem jest zinterpretowanie tego, co reprezentują utworzone grupy i nadanie im nazwy. Interpretacja jest trudnym zadaniem i często wyniki są niejasne. W tym celu konieczne jest przeanalizowanie statystyk opisowych klastrów oraz wykorzystanie różnych wykresów i kompozycji map. Bardzo przydatna jest także znajomość właściwości spektralnych obiektów.
Spróbujmy zatem zinterpretować uzyskane skupienia, korzystając z wykresu pudełkowego. Największe możliwości wizualizacji danych dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”. Wymieniony pakiet wymaga przygotowania zbioru danych do odpowiedniej postaci, tj. dane muszą być przedstawione jako ramka danych w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej konwersji można dokonać w prosty sposób używając pakietu tidyr.
library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danychJak zauważyliśmy wcześniej, nasz zbiór danych jest dość duży i nie ma potrzeby prezentowania wszystkich danych. Możemy to zrobić efektywniej, używając wcześniej wylosowanej próby. Połączmy zatem wylosowane wiersze z macierzy z odpowiadającymi im klastrami (cbind()). Następnie macierz zamienimy na ramkę danych (as.data.frame()).
stats = cbind(mat_omit[idx, ], klaster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats) blue green red nir klaster
1 0.0728 0.0940 0.1100 0.2035 1
2 0.0121 0.0302 0.0217 0.3364 3
3 0.0258 0.0606 0.0341 0.3616 3
4 0.0166 0.0406 0.0253 0.3157 3
5 0.0770 0.1080 0.2156 0.2633 1
6 0.0368 0.0670 0.0572 0.3219 3Wyświetlone dane mają formę szeroką (każdy kanał spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu wykorzystamy funkcję pivot_longer().
stats = pivot_longer(stats, cols = 1:4, names_to = "kanal", values_to = "wartosc")Dla formalności możemy jeszcze zmienić typ danych (klastrów i kanałów) na kategoryczny (factor). W praktyce związane jest to z uproszczeniem struktury danych (przejście ze skali ilorazowej do nominalnej).
stats$klaster = factor(stats$klaster)
stats$kanal = factor(stats$kanal)
head(stats)# A tibble: 6 × 3
klaster kanal wartosc
<fct> <fct> <dbl>
1 1 blue 0.0728
2 1 green 0.094
3 1 red 0.11
4 1 nir 0.204
5 3 blue 0.0121
6 3 green 0.0302Ramka danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.
ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
geom_boxplot()
Zmieńmy kilka domyślnych parametrów żeby poprawić odbiór ryciny.
etykiety = c("Niebieski", "Zielony", "Czerwony", "Bliska\npodczerwień")
ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
geom_boxplot(show.legend = FALSE) +
scale_x_discrete(limits = kanaly, labels = etykiety) +
scale_fill_manual(values = kolory) +
facet_wrap(vars(klaster)) +
xlab("Kanał") +
ylab("Odbicie") +
theme_light()
Na podstawie powyższego wykresu możemy przeanalizować właściwości spektralne klastrów, a tym samym zinterpretować, jakie obiekty reprezentują.
Finalna mapa
Ostatnim etapem jest stworzenie mapy klasyfikacyjnej pokrycia terenu na podstawie otrzymanego wektora z klastrami (mdl$cluster). Na początku musimy przygotować pusty wektor składający się z całkowitej liczby komórek rastra. Można to sprawdzić za pomocą funkcji ncell(). W naszym przypadku jest to 6 milionów komórek.
# przygotuj pusty wektor
wek = rep(NA, ncell(r))Następnie musimy przypisać nasze grupy w wektorze w odpowiednie miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych wartości można odwołać się przez funkcję complete.cases().
# zastąp tylko te wartości, które nie są NA
wek[complete.cases(mat)] = mdl$cluster W ostatnim kroku należy skopiować metadane obiektu r, ale tylko z jedną warstwą, i przypisać mu wartości wektora wek.
# stwórz nowy raster
clustering = rast(r, nlyrs = 1, vals = wek)Zaprezentujmy wynik grupowania na mapie, używając odpowiednich kolorów i nazw klastrów.
kolory = c("#d9d9d9", "#086209", "#2fbd2f")
kategorie = c("odkryta gleba", "lasy/woda", "roślinność")
plot(clustering, col = kolory, type = "classes", levels = kategorie,
mar = c(3, 3, 3, 7))
Jeśli wynik jest zadowalający, możemy zapisać go za pomocą funkcji writeRaster(). Taki plik można później wczytać w R lub innym programie obsługującym dane przestrzenne (np. QGIS). Dodatkowo, w przypadku danych kategorycznych, podczas zapisu warto ustawić typ danych jako Byte / INT1U (pod warunkiem, że liczba kategorii nie przekracza 255).
writeRaster(clustering, "clustering.tif", datatype = "INT1U")Zadanie
7. Pobierz scenę satelitarną o niskim zachmurzeniu i dotnij jej zasięg do wybranego powiatu. Wykonaj klasteryzację metodą kmeans oraz inną wybraną. Następnie dokonaj walidacji otrzymanych wyników wykorzystując wskaźnik silhouette. Przygotuj również wykres pudełkowy przedstawiający zmienność powstałych klastrów. Finalnie, zaprezentuj wynik klasteryzacji na mapie (dobierz odpowiedni schemat kolorów) oraz kompozycję RGB.